{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to atomman: Basic support and analysis tools\n",
    "\n",
    "__Lucas M. Hale__, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), _Materials Science and Engineering Division, NIST_.\n",
    "    \n",
    "[Disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction<a id='section1'></a>\n",
    "\n",
    "This Notebook outlines some of the other tools in atomman that provide basic support features and simple analysis of the atomistic systems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Library Imports**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "atomman version = 1.4.0\n",
      "Notebook executed on 2021-08-05\n"
     ]
    }
   ],
   "source": [
    "# Standard Python libraries\n",
    "import os\n",
    "from io import open\n",
    "from copy import deepcopy\n",
    "import datetime\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np\n",
    "\n",
    "# https://pandas.pydata.org/\n",
    "import pandas as pd\n",
    "\n",
    "# https://github.com/usnistgov/atomman\n",
    "import atomman as am\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# Show atomman version\n",
    "print('atomman version =', am.__version__)\n",
    "\n",
    "# Show date of Notebook execution\n",
    "print('Notebook executed on', datetime.date.today())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct a demonstration 2x2x2 diamond cubic silicon system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64\n"
     ]
    }
   ],
   "source": [
    "a = uc.set_in_units(5.431, 'angstrom')\n",
    "ucell = am.load('prototype', 'A4--C--dc', a=a, symbols='Si')\n",
    "\n",
    "system = ucell.supersize(2,2,2)\n",
    "\n",
    "print(system.natoms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Elastic constants <a id='section2'></a>\n",
    "\n",
    "The full elastic constants tensor for a given crystal can be represented with the atomman.ElasticConstants class.  The values in an ElasticConstants object can be set and retrieved in a variety of formats and transformed to other Cartesian coordinate systems. \n",
    "\n",
    "See the [3.1. ElasticConstants class Jupyter Notebook](3.1._ElasticConstants_class.html) for more details and a full description of all of the class methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define an ElasticConstants object for diamond cubic silicon\n",
    "# values taken from http://www.ioffe.ru/SVA/NSM/Semicond/Si/mechanic.html\n",
    "C11 = uc.set_in_units(16.60 * 10**11, 'dyn/cm^2')\n",
    "C12 = uc.set_in_units( 6.40 * 10**11, 'dyn/cm^2')\n",
    "C44 = uc.set_in_units( 7.96 * 10**11, 'dyn/cm^2')\n",
    "\n",
    "C = am.ElasticConstants(C11=C11, C12=C12, C44=C44)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cij (GPa) =\n",
      "[[166.   64.   64.    0.    0.    0. ]\n",
      " [ 64.  166.   64.    0.    0.    0. ]\n",
      " [ 64.   64.  166.    0.    0.    0. ]\n",
      " [  0.    0.    0.   79.6   0.    0. ]\n",
      " [  0.    0.    0.    0.   79.6   0. ]\n",
      " [  0.    0.    0.    0.    0.   79.6]]\n"
     ]
    }
   ],
   "source": [
    "# Get 6x6 Cij Voigt representation of elastic constants in GPa\n",
    "print('Cij (GPa) =')\n",
    "print(uc.get_in_units(C.Cij, 'GPa'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Relative distances between atoms<a id='section3'></a>\n",
    "\n",
    "There are a few built-in tools for investigating the relative positions between atoms of the same and different systems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1. System.dvect()\n",
    "\n",
    "The System.dvect() method computes the shortest vector(s) between two points or list of points within the System taking into account the System's periodic dimensions.\n",
    "\n",
    "Parameters\n",
    "\n",
    "- __pos_0__ (*numpy.ndarray or index*) Absolute Cartesian vector position(s) to use as reference point(s). If the value can be used as an index, then self.atoms.pos[pos_0] is taken.\n",
    "- __pos_1__ (*numpy.ndarray or index*) Absolute Cartesian vector position(s) to find relative to pos_0.  If the value can be used as an index, then self.atoms.pos[pos_1] is taken."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 4.07325  4.07325 -4.07325]\n"
     ]
    }
   ],
   "source": [
    "# Calculate shortest vector between atoms 1 and 60\n",
    "print(system.dvect(1, 60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2. System.dmag()\n",
    "\n",
    "The System dmag() method computes the magnitude of the shortest vector(s) between two points or list of points within the System taking into account the System's periodic dimensions.  This is identical to computing dvect above, then finding the magnitude of those vectors, but should be faster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[8.66025404 5.95297241 5.95297241 5.95297241 6.30856205 3.87088054\n",
      " 3.87088054 3.87088054 7.08419092 6.33398788 6.33398788 3.25939281\n",
      " 5.45266877 2.2175116  5.86626957 5.86626957 7.08419092 6.33398788\n",
      " 3.25939281 6.33398788 5.45266877 5.86626957 2.2175116  5.86626957\n",
      " 5.03701519 6.69334927 3.91218142 3.91218142 4.43455051 4.93424363\n",
      " 4.93424363 7.33774633 7.08419092 3.25939281 6.33398788 6.33398788\n",
      " 5.45266877 5.86626957 5.86626957 2.2175116  5.03701519 3.91218142\n",
      " 6.69334927 3.91218142 4.43455051 4.93424363 7.33774633 4.93424363\n",
      " 5.03701519 3.91218142 3.91218142 6.69334927 4.43455051 7.33774633\n",
      " 4.93424363 4.93424363 0.7465139  4.4706471  4.4706471  4.4706471\n",
      " 3.09820588 6.6163557  6.6163557  6.6163557 ]\n"
     ]
    }
   ],
   "source": [
    "# Calculate shortest distance between position [5., 5., 5.] and all atoms in system\n",
    "dmags = system.dmag([5.0, 5.0, 5.0], range(system.natoms))\n",
    "print(dmags)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3. displacement()\n",
    "\n",
    "The atomman.displacement() function compares two systems with the same number of atoms and calculates the vector differences between all atoms with the same atomic id's. The vectors returned are the shortest vectors after taking periodic boundaries in consideration, i.e. it uses dvect().\n",
    "\n",
    "Parameters\n",
    "\n",
    "- __system_0__ (*atomman.System*) The initial system to calculate displacements from.\n",
    "- __system_1__ (*atomman.System*) The final system to calculate displacements to.\n",
    "- __box_reference__ (*str or None*) Specifies which system's boundary conditions to use.\n",
    "    - 'initial' uses system_0's box and pbc.\n",
    "    - 'final' uses system_1's box and pbc (Default).\n",
    "    - None computes the straight difference between the positions without accounting for periodic boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.99603857 1.43863724 0.591721  ]\n",
      " [1.72792262 0.61085361 1.70718691]\n",
      " [0.15549565 0.92514288 2.61097281]\n",
      " [1.41569634 0.90522856 1.68847287]\n",
      " [1.81424267 0.56808247 2.18561857]\n",
      " [2.72553133 2.13598235 2.74308199]\n",
      " [1.29520816 2.46282183 2.21613065]\n",
      " [2.90929418 1.73635724 0.41017354]\n",
      " [1.76887711 2.85857528 2.65185785]\n",
      " [0.26059387 1.02830317 0.77041408]\n",
      " [1.27722693 2.63985422 1.0868175 ]\n",
      " [2.64540753 2.67385585 0.48848634]\n",
      " [2.88070417 2.91077294 1.92749443]\n",
      " [0.69530365 2.35475947 0.46623167]\n",
      " [1.84750822 2.9051281  0.3921404 ]\n",
      " [0.64373951 0.24626612 0.8194635 ]\n",
      " [1.89782266 2.66502855 0.46479666]\n",
      " [1.7791236  0.22940824 0.36557604]\n",
      " [2.05966212 0.88053363 0.27463514]\n",
      " [0.99311977 0.48984318 0.38610387]\n",
      " [0.50177108 1.85573905 2.17549783]\n",
      " [2.00151313 1.11888152 1.91759207]\n",
      " [2.82831279 0.51041983 0.43659667]\n",
      " [0.16477686 1.31245116 1.84300738]\n",
      " [1.14831696 0.89599073 0.94318064]\n",
      " [0.40920217 2.69268077 0.92314081]\n",
      " [1.92433531 0.43145832 2.57552459]\n",
      " [1.29957253 0.52358684 2.54466927]\n",
      " [0.37003548 0.80561067 2.4623484 ]\n",
      " [0.71855465 2.70688489 2.43903939]\n",
      " [0.04998238 0.99324663 0.62647591]\n",
      " [1.40175192 1.91228962 0.14133925]\n",
      " [0.873906   0.73383899 1.67054587]\n",
      " [2.91306592 2.15287323 2.97282268]\n",
      " [2.48032709 1.18365632 0.97460698]\n",
      " [2.87915574 1.52372723 0.02143964]\n",
      " [0.18544469 0.21583727 2.94958799]\n",
      " [2.9889433  0.34263071 2.15844326]\n",
      " [2.68588446 1.59649889 1.6853371 ]\n",
      " [2.85236091 1.28573587 2.12574881]\n",
      " [0.70070409 0.79190137 0.71048702]\n",
      " [1.50671995 2.64742722 1.4868234 ]\n",
      " [0.96455881 0.20132281 2.3701744 ]\n",
      " [2.44793886 0.87574734 1.39985676]\n",
      " [1.45910007 1.68904979 0.92585075]\n",
      " [0.16827792 0.25317677 2.4503324 ]\n",
      " [2.04055577 0.91133765 1.99867685]\n",
      " [0.13699851 0.32462108 0.5213217 ]\n",
      " [0.8698135  2.96181797 0.83466607]\n",
      " [0.70480033 0.79672071 1.19483379]\n",
      " [2.99017617 0.93501253 0.74366069]\n",
      " [1.03433172 0.09130963 0.38036951]\n",
      " [0.19494491 0.61943101 2.85282479]\n",
      " [1.3757641  1.625531   0.89416238]\n",
      " [0.19562443 0.77011829 1.51996765]\n",
      " [2.2330411  0.94005377 2.51933319]\n",
      " [1.29213491 0.70965131 0.31615813]\n",
      " [2.77540236 0.85584758 0.64817984]\n",
      " [0.20847362 2.39985079 1.73414601]\n",
      " [0.78406961 2.0580142  1.25450911]\n",
      " [1.98109968 1.63475414 0.67012227]\n",
      " [0.37131018 0.41426808 2.09247931]\n",
      " [0.83840685 0.07476563 0.8252753 ]\n",
      " [1.06777949 2.27272563 0.38276492]]\n"
     ]
    }
   ],
   "source": [
    "# Copy system and randomly displace atoms\n",
    "system2 = deepcopy(system)\n",
    "system2.atoms.pos += 3 * np.random.rand(system.natoms, 3)\n",
    "system2.wrap()\n",
    "\n",
    "# Show displacement between the two systems\n",
    "print(am.displacement(system, system2))       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4. System.neighborlist()\n",
    "\n",
    "A list of neighbor atoms within a cutoff can be constructed using the System.neighborlist() method.  The list of neighbors is returned as an atomman.NeighborList object.\n",
    "\n",
    "See the [3.2. NeighborList class Jupyter Notebook](3.2._NeighborList_class.html) for more details on how the list is calculated and can be used.\n",
    "\n",
    "Parameters\n",
    "        \n",
    "- __cutoff__ (*float, optional*) Radial cutoff distance for identifying neighbors.  Must be given if model is not given.\n",
    "- __model__ (*str or file-like object, optional*) Gives the file path or content to load.  If given, no other parameters are allowed.\n",
    "- __initialsize__ (*int, optional*) The number of neighbor positions to initially assign to each atom.  Default value is 20.\n",
    "- __deltasize__ (*int, optional*) Specifies the number of extra neighbor positions to allow each atom when the number of neighbors exceeds the underlying array size.  Default value is 10.\n",
    "            \n",
    "Returns\n",
    "        \n",
    "- (*atomman.NeighborList*) The compiled list of neighbors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify neighbors within 3 angstroms\n",
    "neighbors = system.neighborlist(cutoff=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The coordinataion numbers for the atoms can be retrieved with coord."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average coordination = 4.0\n"
     ]
    }
   ],
   "source": [
    "# Show average atomic coordination\n",
    "print('Average coordination =', neighbors.coord.mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Index selection on the NeighborList object will return the indices of the neighbor atoms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Neighbors of atom 6 = [ 2 11 33 40]\n"
     ]
    }
   ],
   "source": [
    "# List neighbor atoms of atom 6\n",
    "print('Neighbors of atom 6 =', neighbors[6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The neighbor indices can then be used to filter other properties to only focus on an atom's neighbors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.35169198 2.35169198 2.35169198 2.35169198]\n"
     ]
    }
   ],
   "source": [
    "# List the dmag values between atom 9 and its neighbors\n",
    "print(system.dmag(9, neighbors[9]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Region selectors<a id='section4'></a>\n",
    "\n",
    "*Added version 1.3.0*\n",
    "\n",
    "A number of geometric shape definitions are available in the atomman.region submodule to help identify regions in space above/below planes or inside/outside of regions. These are useful for constructing systems by slicing away atoms to create nanostructures, or for performing analysis on only select regions.\n",
    "\n",
    "See the [3.3. Region selectors Jupyter Notebook](3.3._Region_selectors.html) for more details and a list of all available shapes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64 atoms in system\n",
      "24 atoms above plane\n",
      "40 atoms below plane\n",
      "11 atoms inside sphere\n"
     ]
    }
   ],
   "source": [
    "# Define a plane normal to the y axis and positioned halfway across system\n",
    "plane = am.region.Plane([0,1,0], system.box.bvect / 2)\n",
    "\n",
    "# Count number of atoms in system, and above/below plane\n",
    "print(f'{system.natoms} atoms in system')\n",
    "\n",
    "abovecount = np.sum(plane.above(system.atoms.pos))\n",
    "print(f'{abovecount} atoms above plane')\n",
    "\n",
    "belowcount = np.sum(plane.below(system.atoms.pos))\n",
    "print(f'{belowcount} atoms below plane')\n",
    "\n",
    "# Define a sphere centered at [0,0,0] with radius = 6\n",
    "sphere = am.region.Sphere([0,0,0], 6)\n",
    "\n",
    "# Count atoms inside sphere\n",
    "insidecount = np.sum(sphere.inside(system.atoms.pos))\n",
    "print(f'{insidecount} atoms inside sphere')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Basic tools<a id='section5'></a>\n",
    "\n",
    "This lists some of the other basic tools and features in atomman."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1. Atomic information\n",
    "\n",
    "- __atomman.tools.atomic_number()__ returns the atomic number associated with an element's atomic symbol.  \n",
    "- __atomman.tools.atomic_symbol()__ returns the elemental symbol associated with an given atomic number.\n",
    "- __atomman.tools.atomic_mass()__ returns the atomic mass of an element or isotope. The atom can be identified with atomic number or atomic/isotope symbol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26\n",
      "Fe\n",
      "55.845\n",
      "55.845\n"
     ]
    }
   ],
   "source": [
    "# Get atomic number for an atomic symbol\n",
    "num = am.tools.atomic_number('Fe')\n",
    "print(num)\n",
    "\n",
    "# Get atomic symbol for an atomic number\n",
    "symbol = am.tools.atomic_symbol(num)\n",
    "print(symbol)\n",
    "\n",
    "# Get atomic mass for an atomic symbol\n",
    "mass = am.tools.atomic_mass(symbol)\n",
    "print(mass)\n",
    "\n",
    "# Get atomic mass for an atomic number\n",
    "mass = am.tools.atomic_mass(num)\n",
    "print(mass)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "25.986891904\n"
     ]
    }
   ],
   "source": [
    "# Get atomic mass for an isotope\n",
    "mass = am.tools.atomic_mass('Al-26')\n",
    "print(mass)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.2. axes_check()\n",
    "\n",
    "The axes_check() function is useful when working in Cartesian systems. Given a (3,3) array representing three 3D Cartesian vectors:\n",
    "\n",
    "- The three vectors are checked that they are orthogonal and right-handed.\n",
    "- The corresponding array of unit vectors are returned. This can then be used for crystal transformations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.70710678  0.          0.70710678]\n",
      " [ 0.70710678  0.          0.70710678]\n",
      " [ 0.          1.          0.        ]]\n"
     ]
    }
   ],
   "source": [
    "axes = [[-1, 0, 1],\n",
    "        [ 1, 0, 1],\n",
    "        [ 0, 1, 0]]\n",
    "print(am.tools.axes_check(axes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.3. filltemplate()\n",
    "\n",
    "The filltemplate() function takes a template and fills in values for delimited template variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "My friend Charlie really likes to use templates to program, says that they are delicious!\n"
     ]
    }
   ],
   "source": [
    "madlibs = \"My friend <name> really likes to use templates to <verb>, says that they are <adjective>!\"\n",
    "s_delimiter = '<'\n",
    "e_delimiter = '>'\n",
    "\n",
    "terms = {}\n",
    "terms['name'] = 'Charlie'\n",
    "terms['verb'] = 'program'\n",
    "terms['adjective'] = 'delicious'\n",
    "\n",
    "print(am.tools.filltemplate(madlibs, terms, s_delimiter, e_delimiter))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.4. indexstr()\n",
    "\n",
    "Iterates through all indicies of an array with a given shape, returning both the numeric index and a string representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "index -> (0, 0) , istr -> '[0][0]'\n",
      "index -> (0, 1) , istr -> '[0][1]'\n",
      "index -> (1, 0) , istr -> '[1][0]'\n",
      "index -> (1, 1) , istr -> '[1][1]'\n",
      "index -> (2, 0) , istr -> '[2][0]'\n",
      "index -> (2, 1) , istr -> '[2][1]'\n"
     ]
    }
   ],
   "source": [
    "for index, istr in am.tools.indexstr((3,2)):\n",
    "    print('index ->', repr(index), ', istr ->', repr(istr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.5. uber_open_rmode\n",
    "\n",
    "uber_open_rmode is a context manager that allows for similar reading of content from a file or from a string variable. It equivalently handles:\n",
    "    \n",
    "- str path name to a file\n",
    "- str content\n",
    "- open file-like object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b'Here I am, read me!'\n",
      "b'Here I am, read me!'\n",
      "b'Here I am, read me!'\n"
     ]
    }
   ],
   "source": [
    "# Define str and save to file\n",
    "text = 'Here I am, read me!'\n",
    "fname = 'text.txt'\n",
    "with open(fname, 'w') as f:\n",
    "    f.write(text)\n",
    "\n",
    "# Use uber_open_rmode on text\n",
    "with am.tools.uber_open_rmode(text) as f:\n",
    "    print(f.read())\n",
    "    \n",
    "# Use uber_open_rmode on file path\n",
    "with am.tools.uber_open_rmode(fname) as f:\n",
    "    print(f.read())\n",
    "    \n",
    "# Use uber_open_rmode on file-like object\n",
    "with open(fname, 'rb') as fobject:\n",
    "    with am.tools.uber_open_rmode(fobject) as f:\n",
    "        print(f.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.6. vect_angle()\n",
    "\n",
    "The vect_angle() function returns the angle between two vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Angle between [-0.34828678 -0.74741765  0.60540275] and [-0.67224036  0.72470423 -0.14586779] =\n",
      "112.78423936770139 degrees\n",
      "1.9684563213237967 radians\n"
     ]
    }
   ],
   "source": [
    "vect1 = 2*np.random.rand(3)-1\n",
    "vect2 = 2*np.random.rand(3)-1\n",
    "\n",
    "print('Angle between', vect1, 'and', vect2, '=')\n",
    "print(am.tools.vect_angle(vect1, vect2), 'degrees')\n",
    "print(am.tools.vect_angle(vect1, vect2, 'radian'), 'radians')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.7. duplicates_allclose()\n",
    "\n",
    "Determine duplicates in dataframe based on tolerances.  The implementation\n",
    "first uses pandas.DataFrame.duplicated on the `dcols` argument with\n",
    "`keep=False` to keep all duplicates.  The duplicate sub-dataframe is then\n",
    "sorted on both `dcols` and `fcols`.  A diff between each row is then done\n",
    "on the sorted duplicates dataframe.  The float values are then checked for\n",
    "their tolerances.\n",
    "\n",
    "Note: False duplicates may be identified if tolerance ranges overlap.\n",
    "Consider dataframe with rows 1,2,3.  If row 2 matches row 1 within the\n",
    "tolerances, and row 3 matches row 2 within the tolerances, both rows 2 and\n",
    "3 will be labeled as tolerances even if row 3 does not match row 1 within\n",
    "the tolerances.\n",
    "\n",
    "Parameters\n",
    "- __dataframe__ (*pandas.DataFrame*) The dataframe to search for duplicates\n",
    "- __dcols__ (*list*) The column names that are tested for exact duplicates.\n",
    "- __fcols__ (*dict*) The column names (keys) that are tested using absolute tolerances (values).\n",
    "\n",
    "Returns\n",
    "- (*list of bool of length nrows*) False for first occurrence of checked values, True for subsequent duplicates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>0</td>\n",
       "      <td>1.00001</td>\n",
       "      <td>Same</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>1.00002</td>\n",
       "      <td>Diff</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>3.00000</td>\n",
       "      <td>Same</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>1.00000</td>\n",
       "      <td>Same</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         A     B\n",
       "0  1.00001  Same\n",
       "1  1.00002  Diff\n",
       "2  3.00000  Same\n",
       "3  1.00000  Same"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Generate test DataFrame\n",
    "df = pd.DataFrame({'A':[1.00001, 1.00002, 3.0000, 1.000000], 'B':['Same', 'Diff', 'Same', 'Same']})\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <td>1</td>\n",
       "      <td>1.00002</td>\n",
       "      <td>Diff</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>2</td>\n",
       "      <td>3.00000</td>\n",
       "      <td>Same</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <td>3</td>\n",
       "      <td>1.00000</td>\n",
       "      <td>Same</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         A     B\n",
       "1  1.00002  Diff\n",
       "2  3.00000  Same\n",
       "3  1.00000  Same"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Show unique values \n",
    "df[~am.tools.duplicates_allclose(df, dcols=['B'], fcols={'A':1e-4})]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.8. Miller index conversions\n",
    "\n",
    "- __atomman.tools.miller.vector3to4(indices)__ converts vectors from three-term Miller indices to four-term Miller-Bravais indices for hexagonal systems. *Updated version 1.2.6*: now returns vectors with magnitudes consistent withe given vectors rather than rescaling to the smallest integer representations.\n",
    "- __atomman.tools.miller.vector4to3(indices)__ converts vectors from four-term Miller-Bravais indices to three-term Miller indices. *Updated version 1.2.6*: now returns vectors with magnitudes consistent withe given vectors rather than rescaling to the smallest integer representations.\n",
    "- __atomman.tools.miller.plane3to4(indices)__ converts planes from three-term Miller indices to four-term Miller-Bravais indices for hexagonal systems. *Added version 1.2.8*\n",
    "- __atomman.tools.miller.plane4to3(indices)__ converts planes from four-term Miller-Bravais indices to three-term Miller indices. *Added version 1.2.8*\n",
    "- __atomman.tools.miller.vector_crystal_to_cartesian(indices, box)__ converts Miller and Miller-Bravais indices to Cartesian vectors based on a supplied box. *Updated version 1.2.6* renamed from vectortocartesian for consistency.\n",
    "- __atomman.tools.miller.plane_crystal_to_cartesian(indices, box)__ converts Miller and Miller-Bravais plane indices to Cartesian normal vectors based on a supplied box. The method uses the definition of the crystal planes to identify two in-plane crystal vectors, converts them to Cartesian, and obtains the plane normal as the cross product of the two vectors. *Added version 1.3.2* \n",
    "- __atomman.tools.miller.vector_primitive_to_conventional(indices, setting)__ converts vectors relative to a primitive unit cell to a conventional unit cell in the given setting (p, a, b, c, i or f). *Added version 1.2.6*\n",
    "- __atomman.tools.miller.vector_conventional_to_primitive(indices, setting)__ converts vectors relative to a conventional unit cell in the given setting (p, a, b, c, i or f) to a primitive unit cell. *Added version 1.2.6*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.  1. -2.  3.]\n",
      "[3. 3. 0.]\n"
     ]
    }
   ],
   "source": [
    "# Test single vector conversions\n",
    "print(am.tools.miller.vector3to4(np.array([3,3,3])))\n",
    "print(am.tools.miller.vector4to3(np.array([1,1,-2,0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1  0  2]\n",
      " [ 3 -4  4]\n",
      " [ 1 -2  1]]\n",
      "\n",
      "[[-0.66666667  0.33333333  0.33333333  2.        ]\n",
      " [ 3.33333333 -3.66666667  0.33333333  4.        ]\n",
      " [ 1.33333333 -1.66666667  0.33333333  1.        ]]\n",
      "\n",
      "[[-1.  0.  2.]\n",
      " [ 3. -4.  4.]\n",
      " [ 1. -2.  1.]]\n"
     ]
    }
   ],
   "source": [
    "# Generate random uvw crystal indices\n",
    "indices = np.random.randint(-5,6, (3,3))\n",
    "print(indices)\n",
    "print()\n",
    "\n",
    "# Convert to hexagonal uvtw's\n",
    "indices = am.tools.miller.vector3to4(indices)\n",
    "print(indices)\n",
    "print()\n",
    "\n",
    "# Convert back to uvw's and see that values are recovered\n",
    "indices = am.tools.miller.vector4to3(indices)\n",
    "print(indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.  3. -6.  3.]\n",
      "[1. 1. 0.]\n"
     ]
    }
   ],
   "source": [
    "# Test single plane conversions\n",
    "print(am.tools.miller.plane3to4(np.array([3,3,3])))\n",
    "print(am.tools.miller.plane4to3(np.array([1,1,-2,0])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-2  4  0]\n",
      " [-1 -3  5]\n",
      " [ 3  3 -5]]\n",
      "\n",
      "[[-2.  4. -2.  0.]\n",
      " [-1. -3.  4.  5.]\n",
      " [ 3.  3. -6. -5.]]\n",
      "\n",
      "[[-2.  4.  0.]\n",
      " [-1. -3.  5.]\n",
      " [ 3.  3. -5.]]\n"
     ]
    }
   ],
   "source": [
    "# Generate random hkl crystal indices\n",
    "indices = np.random.randint(-5,6, (3,3))\n",
    "print(indices)\n",
    "print()\n",
    "\n",
    "# Convert to hexagonal hkil's\n",
    "indices = am.tools.miller.plane3to4(indices)\n",
    "print(indices)\n",
    "print()\n",
    "\n",
    "# Convert back to hkl's and see that values are recovered\n",
    "indices = am.tools.miller.plane4to3(indices)\n",
    "print(indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.51        0.          0.        ]\n",
      " [-1.255       2.17372376  0.        ]\n",
      " [ 0.          0.          4.07      ]]\n",
      "\n",
      "[[ 2.51        0.          0.        ]\n",
      " [-1.255       2.17372376  0.        ]\n",
      " [ 0.          0.          4.07      ]]\n"
     ]
    }
   ],
   "source": [
    "# Define a hexagonal box\n",
    "a = uc.set_in_units(2.51, 'angstrom')\n",
    "c = uc.set_in_units(4.07, 'angstrom')\n",
    "box = am.Box(a=a, b=a, c=c, gamma=120)\n",
    "\n",
    "# Pass Miller indices\n",
    "indices = [[1,0,0],\n",
    "           [0,1,0],\n",
    "           [0,0,1]]\n",
    "print(am.tools.miller.vector_crystal_to_cartesian(indices, box))\n",
    "print()\n",
    "\n",
    "# Pass equivalent Miller-Bravais indices\n",
    "indices = [[ 2/3,-1/3,-1/3, 0],\n",
    "           [-1/3, 2/3,-1/3, 0],\n",
    "           [   0,   0,   0, 1]]\n",
    "print(am.tools.miller.vector_crystal_to_cartesian(indices, box))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0. -0.  1.]\n",
      "[ 0.8660254  0.5       -0.       ]\n",
      "\n",
      "[ 0. -0.  1.]\n",
      "[ 0.8660254  0.5       -0.       ]\n"
     ]
    }
   ],
   "source": [
    "# Define a hexagonal box\n",
    "a = uc.set_in_units(2.51, 'angstrom')\n",
    "c = uc.set_in_units(4.07, 'angstrom')\n",
    "box = am.Box(a=a, b=a, c=c, gamma=120)\n",
    "\n",
    "# Pass Miller plane indices\n",
    "indices = [ 0, 0, 1]\n",
    "print(am.tools.miller.plane_crystal_to_cartesian(indices, box))\n",
    "indices = [ 1, 0, 0]\n",
    "print(am.tools.miller.plane_crystal_to_cartesian(indices, box))\n",
    "print()\n",
    "\n",
    "# Pass equivalent Miller-Bravais indices\n",
    "indices = [ 0, 0, 0, 1]\n",
    "print(am.tools.miller.plane_crystal_to_cartesian(indices, box))\n",
    "indices = [ 1, 0,-1, 0]\n",
    "print(am.tools.miller.plane_crystal_to_cartesian(indices, box))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 2.477,  0.000,  0.000]\n",
      "bvect =  [-0.825,  2.335,  0.000]\n",
      "cvect =  [-0.825, -1.167,  2.023]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 1\n",
      "natypes = 1\n",
      "symbols = (None,)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Define a primitive bcc unit cell box\n",
    "a = uc.set_in_units(2.86, 'angstrom')\n",
    "p_box = am.Box.trigonal(a * 3**0.5 / 2, alpha=109.466666667)\n",
    "p_ucell = am.System(box=p_box)\n",
    "print(p_ucell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "primitive uvws:\n",
      "[[ 0. -1. -1.]\n",
      " [ 1.  1.  0.]\n",
      " [ 1.  0.  1.]]\n",
      "conventional uvws:\n",
      "[[1. 0. 0.]\n",
      " [0. 1. 0.]\n",
      " [0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "# Convert conventional box vectors to primitive vectors\n",
    "a_uvw = am.tools.miller.vector_conventional_to_primitive([1, 0, 0], setting='i')\n",
    "b_uvw = am.tools.miller.vector_conventional_to_primitive([0, 1, 0], setting='i')\n",
    "c_uvw = am.tools.miller.vector_conventional_to_primitive([0, 0, 1], setting='i')\n",
    "p_uvws = np.array([a_uvw, b_uvw, c_uvw])\n",
    "print('primitive uvws:')\n",
    "print(p_uvws)\n",
    "\n",
    "# Convert back to conventional just for consistency\n",
    "print('conventional uvws:')\n",
    "print(am.tools.miller.vector_primitive_to_conventional(p_uvws, setting='i'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 2.860,  0.000,  0.000]\n",
      "bvect =  [-0.000,  2.860,  0.000]\n",
      "cvect =  [-0.000,  0.000,  2.860]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 2\n",
      "natypes = 1\n",
      "symbols = (None,)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   1.430 |   1.430 |   1.430\n"
     ]
    }
   ],
   "source": [
    "# rotate system using p_uvws to get conventional unit cell\n",
    "c_ucell = p_ucell.rotate(p_uvws)\n",
    "print(c_ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.9. Crystal lattice identification\n",
    "\n",
    "There are also a few tests for identifying if a supplied box is consistent with a standard representation of a crystal family unit cell.\n",
    "\n",
    "- __atomman.tools.identifyfamily(box)__ returns str crystal family if box corresponds to a standard crystal representation. Otherwise, returns None.\n",
    "- __atomman.tools.iscubic(box))__ returns bool indicating if box is a standard cubic box.\n",
    "- __atomman.tools.ishexagonal(box))__ returns bool indicating if box is a standard hexagonal box.\n",
    "- __atomman.tools.istetragonal(box))__ returns bool indicating if box is a standard tetragonal box.\n",
    "- __atomman.tools.isrhombohedral(box))__ returns bool indicating if box is a standard rhombohedral box.\n",
    "- __atomman.tools.isorthorhombic(box))__ returns bool indicating if box is a standard orthorhombic box.\n",
    "- __atomman.tools.ismonoclinic(box))__ returns bool indicating if box is a standard monoclinic box.\n",
    "- __atomman.tools.istriclinic(box))__ returns bool indicating if box is a standard triclinic box.\n",
    "\n",
    "All of these functions use the following standard representation criteria:\n",
    "\n",
    "- cubic: \n",
    "    - $a = b = c$\n",
    "    - $\\alpha = \\beta = \\gamma = 90$\n",
    "- hexagonal: \n",
    "    - $a = b \\ne c$\n",
    "    - $\\alpha = \\beta = 90$\n",
    "    - $\\gamma = 120$\n",
    "- tetragonal: \n",
    "    - $a = b \\ne c$\n",
    "    - $\\alpha = \\beta = \\gamma = 90$\n",
    "- rhombohedral:\n",
    "    - $a = b = c$\n",
    "    - $\\alpha = \\beta = \\gamma \\ne 90$\n",
    "- orthorhombic: \n",
    "    - $a \\ne b \\ne c$\n",
    "    - $\\alpha = \\beta = \\gamma = 90$\n",
    "- monoclinic: \n",
    "    - $a \\ne b \\ne c$\n",
    "    - $\\alpha = \\gamma = 90$\n",
    "    - $\\beta \\ne 90$\n",
    "- triclinic: \n",
    "    - $a \\ne b \\ne c$\n",
    "    - $\\alpha \\ne \\beta \\ne \\gamma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "identifyfamily = orthorhombic\n",
      "iscubic =        False\n",
      "ishexagonal =    False\n",
      "istetragonal =   False\n",
      "isrhombohedral = False\n",
      "isorthorhombic = True\n",
      "ismonoclinic =   False\n",
      "istriclinic =    False\n"
     ]
    }
   ],
   "source": [
    "# Define an orthogonal box\n",
    "a = uc.set_in_units(2.51, 'angstrom')\n",
    "b = uc.set_in_units(3.13, 'angstrom')\n",
    "c = uc.set_in_units(4.07, 'angstrom')\n",
    "box = am.Box(a=a, b=b, c=c)\n",
    "\n",
    "print('identifyfamily =', am.tools.identifyfamily(box))\n",
    "print('iscubic =       ', am.tools.iscubic(box))\n",
    "print('ishexagonal =   ', am.tools.ishexagonal(box))\n",
    "print('istetragonal =  ', am.tools.istetragonal(box))\n",
    "print('isrhombohedral =', am.tools.isrhombohedral(box))\n",
    "print('isorthorhombic =', am.tools.isorthorhombic(box))\n",
    "print('ismonoclinic =  ', am.tools.ismonoclinic(box))\n",
    "print('istriclinic =   ', am.tools.istriclinic(box))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "identifyfamily = None\n"
     ]
    }
   ],
   "source": [
    "# Define a non-standard tetragonal box with a=c!=b\n",
    "box = am.Box(a=a, b=b, c=a)\n",
    "print('identifyfamily =', am.tools.identifyfamily(box))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.10. compositionstr()\n",
    "\n",
    "*Added version 1.2.7*\n",
    "\n",
    "Takes a list of symbols and the counts for each and returns a reduced composition string.  Used by System.composition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Composition = Al2Si5\n"
     ]
    }
   ],
   "source": [
    "symbols = ['Si', 'Al', 'Si']\n",
    "counts = [500, 1000, 2000]\n",
    "print('Composition =', am.tools.compositionstr(symbols, counts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**File Cleanup**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.remove('text.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
